FUNCTION CONV, x, y
; Compute the convolution of x and y by using FFT.
; x and y are float type arrays.

; Determine the size and type of x and y.
tx = size(x)
ty = size(y)
dx = tx[0] ;Dimension of x
dy = ty[0] ;Dimension of y

; Check dimensions.
IF (dx NE dy) OR (dx GT 2) OR $
   (dy GT 2) THEN BEGIN
    print, "   Dimension ERROR!"
    return, 0
ENDIF

; Check x.
IF (dx EQ 1) AND (tx[2] NE 4) OR $ 
   (dx EQ 2) AND (tx[3] NE 4) THEN BEGIN
    print, "   x's type is not FLOAT. ERROR!"
    return, 0
ENDIF

; Check y.
IF (dy EQ 1) AND (ty[2] NE 4) OR $ 
   (dy EQ 2) AND (ty[3] NE 4) THEN BEGIN
    print, "   y's type is not FLOAT. ERROR!"
    return, 0
ENDIF

; Determine the size of FFT buffer.
lx = tx[1] + ty[1] - 1
lbfx = 2;
WHILE (lbfx LT lx) DO lbfx *= 2
 
IF (dx EQ 2) THEN BEGIN
    ly = tx[2] + ty[2] - 1
    lbfy = 2
    WHILE (lbfy LT ly) DO lbfy *= 2
ENDIF

; Allocate buffer.
IF (dx EQ 1) THEN BEGIN
    bfx = fltarr(lbfx,1)
    bfy = fltarr(lbfx,1)
    
    ; copy the data.
    bfx[0:tx[1]-1] = x
    bfy[0:ty[1]-1] = y
ENDIF ELSE BEGIN
    bfx = fltarr(lbfx, lbfy)
    bfy = fltarr(lbfx, lbfy)

    ; Copy the data.
    bfx[0:tx[1]-1,0:tx[2]-1] = x
    bfy[0:ty[1]-1,0:ty[2]-1] = y
ENDELSE

; Calculate the convolution.
fbfx = fft(bfx)
fbfy = fft(bfy)

fbfxy = fbfx * fbfy
bfxy = abs(fft(fbfxy,/INVERSE))*lbfx

IF (dx EQ 2) THEN bfxy *= lbfy

IF (dx EQ 1) THEN BEGIN
    xy = bfxy[0:lx-1]
ENDIF ELSE BEGIN 
    xy = bfxy[0:lx-1,0:ly-1]
ENDELSE

; Return the results.
return, xy

END
